Análisis de popularidad de canciones en el servicio de streaming Spotify

library(rstan)
library(ggplot2)
library(gridExtra)
library(tidyverse)
library(cmdstanr)
library(rsample)
library(DiagrammeR)
library(kableExtra)

Introducción

De acuerdo con Wikipedia, en la plataforma de Spotify se tiene acceso a más de 100 millones de canciones. Algunas canciones son más populares que otras. Definimos a la variable \(Y\) como un rating de popularidad con \(Y ∈ [0, 100]\). En general, entre más reproducciones recientes tenga una una canción, su rating de popularidad será mayor.

Sabemos que la popularidad de una canción, no siempre es reflejo de su calidad, popularidad a largo plazo, o popularidad afuera de la audiencia de Spotify. Además del artista, hay otras características que hacen popular a una canción, como la “valencia”, que define qué tan alegre suena una pieza.

Para el proyecto nos enfocaremos en entender:

Algunas de las variables que afectan la popularidad de las canciones y en qué grado.
¿Cuál es la popularidad típica de una canción en Spotify?
¿Hasta dónde influye el artista en la popularidad de la canción?
Para un sólo artista, ¿cómo varía la popularidad de sus canciones?
¿Cómo influye la valencia en la popularidad de una canción?

Métodos:

Para nuestro proyecto haremos uso de distintos tipos de modelos:

Para el objetivo inicial, generaremos un modelo lineal generalizado con priors poco informativas haciendo uso de las variables ‘popularity’, ‘duration’, ‘explicit’, ‘danceability’, ‘tempo’ y ‘genre’. Esto nos permitirá conocer los aspectos más relevantes y el impacto de dichas variables en la popularidad de las canciones.
Compararemos el modelo anterior con un modelo con priors informativas.
Modelo con parámetros homogéneos: Modelar la popularidad de las canciones sin tomar en cuenta el artista.
Modelo con parámetros heterogéneos: Modelar la popularidad de las canciones tomando en cuenta a los artistas que las interpretan.
Modelo jerárquico: Una alternativa intermedia donde permitimos que la distribución inicial sobre la popularidad pueda adaptarse a los datos. Y comparar sus distribuciones predictivas con los datos observados, para ver cuál modela mejor y porqué.
Modelo jerárquico con DAG: Haremos un DAG para agregar la variable de valencia al modelo jerárquico.

Datos

Para este proyecto usaremos la base de datos “Spotify Tracks Dataset” de Kaggle, que incluye más de 125 géneros musicales y más de 100 mil canciones. Observamos la estructura de los datos:

datos <- read.csv("datos/datos.csv") 
datos <- as.data.frame(datos)
selected_datos_small <- read.csv("datos/datos_small.csv") 
selected_datos <- read.csv("datos/selected_datos_75.csv") 
glimpse(datos)

## Rows: 114,000
## Columns: 21
## $ id               <int> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,…
## $ track_id         <chr> "5SuOikwiRyPMVoIQDJUgSV", "4qPNDBW1i3p13qLCt0Ki3A", "…
## $ artists          <chr> "Gen Hoshino", "Ben Woodward", "Ingrid Michaelson;ZAY…
## $ album_name       <chr> "Comedy", "Ghost (Acoustic)", "To Begin Again", "Craz…
## $ track_name       <chr> "Comedy", "Ghost - Acoustic", "To Begin Again", "Can'…
## $ popularity       <int> 73, 55, 57, 71, 82, 58, 74, 80, 74, 56, 74, 69, 52, 6…
## $ duration_ms      <int> 230666, 149610, 210826, 201933, 198853, 214240, 22940…
## $ explicit         <chr> "False", "False", "False", "False", "False", "False",…
## $ danceability     <dbl> 0.676, 0.420, 0.438, 0.266, 0.618, 0.688, 0.407, 0.70…
## $ energy           <dbl> 0.4610, 0.1660, 0.3590, 0.0596, 0.4430, 0.4810, 0.147…
## $ key              <int> 1, 1, 0, 0, 2, 6, 2, 11, 0, 1, 8, 4, 7, 3, 2, 4, 2, 1…
## $ loudness         <dbl> -6.746, -17.235, -9.734, -18.515, -9.681, -8.807, -8.…
## $ mode             <int> 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0,…
## $ speechiness      <dbl> 0.1430, 0.0763, 0.0557, 0.0363, 0.0526, 0.1050, 0.035…
## $ acousticness     <dbl> 0.0322, 0.9240, 0.2100, 0.9050, 0.4690, 0.2890, 0.857…
## $ instrumentalness <dbl> 1.01e-06, 5.56e-06, 0.00e+00, 7.07e-05, 0.00e+00, 0.0…
## $ liveness         <dbl> 0.3580, 0.1010, 0.1170, 0.1320, 0.0829, 0.1890, 0.091…
## $ valence          <dbl> 0.7150, 0.2670, 0.1200, 0.1430, 0.1670, 0.6660, 0.076…
## $ tempo            <dbl> 87.917, 77.489, 76.332, 181.740, 119.949, 98.017, 141…
## $ time_signature   <int> 4, 4, 4, 3, 4, 4, 3, 4, 4, 4, 4, 3, 4, 4, 4, 3, 4, 4,…
## $ track_genre      <chr> "acoustic", "acoustic", "acoustic", "acoustic", "acou…

Tenemos 114 mil registros de canciones y 21 características que los representan, entre las cuales destaca el identificador de la canción, el artista, el álbum en el que fue publicada, la popularidad, su duración, entre otras características. Debido al desarrollo del modelado, vamos a ocupar sólo una muestra de 75,000 canciones. Nos enfocaremos en el análisis y modelado de la popularidad, sin embargo, en el anexo 1 de este trabajo se puede consultar la limpieza y adecuación de los datos utilizados.

datos <- read.csv("datos/muestra.csv")
print(summary(datos))

##  artists_numeric       id           track_id           artists         
##  Min.   :   1    Min.   :    25   Length:5016        Length:5016       
##  1st Qu.: 947    1st Qu.: 28783   Class :character   Class :character  
##  Median :1892    Median : 58072   Mode  :character   Mode  :character  
##  Mean   :1879    Mean   : 57685                                        
##  3rd Qu.:2797    3rd Qu.: 86773                                        
##  Max.   :3746    Max.   :113963                                        
##   album_name         track_name          popularity     duration_ms     
##  Length:5016        Length:5016        Min.   : 0.00   Min.   :  30474  
##  Class :character   Class :character   1st Qu.:17.00   1st Qu.: 173870  
##  Mode  :character   Mode  :character   Median :34.00   Median : 213326  
##                                        Mean   :33.07   Mean   : 228133  
##                                        3rd Qu.:50.00   3rd Qu.: 262240  
##                                        Max.   :97.00   Max.   :2851109  
##     explicit        danceability        energy               key        
##  Min.   :0.00000   Min.   :0.0000   Min.   :0.0000195   Min.   : 0.000  
##  1st Qu.:0.00000   1st Qu.:0.4580   1st Qu.:0.4790000   1st Qu.: 2.000  
##  Median :0.00000   Median :0.5800   Median :0.6820000   Median : 6.000  
##  Mean   :0.08174   Mean   :0.5669   Mean   :0.6421708   Mean   : 5.344  
##  3rd Qu.:0.00000   3rd Qu.:0.6940   3rd Qu.:0.8520000   3rd Qu.: 8.000  
##  Max.   :1.00000   Max.   :0.9760   Max.   :1.0000000   Max.   :11.000  
##     loudness            mode       speechiness       acousticness    
##  Min.   :-40.843   Min.   :0.00   Min.   :0.00000   Min.   :0.00000  
##  1st Qu.:-10.142   1st Qu.:0.00   1st Qu.:0.03570   1st Qu.:0.01818  
##  Median : -7.037   Median :1.00   Median :0.04850   Median :0.16450  
##  Mean   : -8.313   Mean   :0.64   Mean   :0.08347   Mean   :0.31228  
##  3rd Qu.: -5.033   3rd Qu.:1.00   3rd Qu.:0.08500   3rd Qu.:0.59300  
##  Max.   :  1.023   Max.   :1.00   Max.   :0.95600   Max.   :0.99600  
##  instrumentalness      liveness         valence           tempo       
##  Min.   :0.00e+00   Min.   :0.0000   Min.   :0.0000   Min.   :  0.00  
##  1st Qu.:0.00e+00   1st Qu.:0.0986   1st Qu.:0.2590   1st Qu.: 99.99  
##  Median :4.29e-05   Median :0.1320   Median :0.4640   Median :121.99  
##  Mean   :1.59e-01   Mean   :0.2134   Mean   :0.4743   Mean   :122.20  
##  3rd Qu.:5.94e-02   3rd Qu.:0.2720   3rd Qu.:0.6873   3rd Qu.:140.11  
##  Max.   :1.00e+00   Max.   :0.9900   Max.   :0.9870   Max.   :213.93  
##  time_signature  track_genre        genre_factor       genre_numeric   
##  Min.   :0.000   Length:5016        Length:5016        Min.   :  1.00  
##  1st Qu.:4.000   Class :character   Class :character   1st Qu.: 29.00  
##  Median :4.000   Mode  :character   Mode  :character   Median : 59.00  
##  Mean   :3.898                                         Mean   : 58.18  
##  3rd Qu.:4.000                                         3rd Qu.: 87.00  
##  Max.   :5.000                                         Max.   :114.00  
##  artists_factor     popularity_artist
##  Length:5016        Min.   : 0.00    
##  Class :character   1st Qu.:18.50    
##  Mode  :character   Median :33.67    
##                     Mean   :33.07    
##                     3rd Qu.:48.00    
##                     Max.   :92.33

Modelos

1) Modelo lineal generalizado

Generamos un modelo lineal con priors poco informativas.

El modelo lineal simple lo expresamos como:

\(y_i = \beta * x_i + \sigma_i\)

Donde:

\(y_i\) es la popularidad de la observación \(i\).
\(x_i\) representa una variable predictora (por ejemplo, ‘duration’) para la observación i.
\(\beta\) representa el coeficiente del modelo.
\(\sigma_i\) es el error aleatorio asociado con la observación \(i\), que asumiremos normalmente distribuido con media cero y varianza constante.

Para nuestro ejercicio requerimos 5 betas para cada una de las variables predictoras que usaremos.

stan_code_lineal <- "
data {
  int<lower=0> N;  
  int<lower=0, upper=100> popularity[N];    
  real<lower=0> duration[N];      
  int<lower=0, upper=1> explicito[N];
  real<lower=0, upper=1> dance[N];  
  real<lower=0, upper=244> tempo[N];     
  int<lower=1, upper=114> genre[N];     
}
parameters {
  real mu_popularity; 
  real<lower=0> sigma_popularity;                            
  real beta_duration;
  real beta_explicito;
  real beta_dance;
  real beta_tempo;
  real beta_genre;                             
}
model {
  // Priors
  mu_popularity ~ normal(0,1);        
  sigma_popularity ~ normal(0,1);         
  beta_duration ~ normal(0,1);        
  beta_explicito ~ normal(0,1);        
  beta_dance ~ normal(0,1);   
  beta_tempo ~ normal(0,1);        
  beta_genre ~ normal(0,1);        
     

  // Likelihood
  for (i in 1:N) {
    popularity[i] ~ normal(mu_popularity + 
                            beta_duration * duration[i] +
                            beta_explicito * explicito[i] +
                            beta_dance * dance[i] +
                            beta_tempo * tempo[i] +
                            beta_genre * genre[i] 
                            , sigma_popularity); 
  }
}
"

# Compilar el modelo
stan_model_lineal <- stan_model(model_code = stan_code_lineal)

# Ajustar el modelo a los datos
fit_lineal <- sampling(stan_model_lineal, 
                       data = list(N = nrow(selected_datos_small), 
                                   popularity = selected_datos_small$popularity,
                                   duration = selected_datos_small$duration_ms,
                                   explicito = selected_datos_small$explicit,
                                   dance = selected_datos_small$danceability,
                                   tempo = selected_datos_small$tempo,
                                   genre = selected_datos_small$genre_numeric), 
                       iter=10000, warmup=3000, seed=123)

## Warning: There were 4676 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded

## Warning: Examine the pairs() plot to diagnose sampling problems

Observamos los resultados

print(fit_lineal)

## Inference for Stan model: anon_model.
## 4 chains, each with iter=10000; warmup=3000; thin=1; 
## post-warmup draws per chain=7000, total post-warmup draws=28000.
## 
##                     mean se_mean   sd    2.5%     25%     50%     75%   97.5%
## mu_popularity       0.29    0.01 1.00   -1.67   -0.39    0.29    0.96    2.27
## sigma_popularity   12.58    0.00 0.47   11.69   12.27   12.57   12.89   13.52
## beta_duration       0.00    0.00 0.00    0.00    0.00    0.00    0.00    0.00
## beta_explicito      1.01    0.01 0.99   -0.92    0.34    1.02    1.68    2.93
## beta_dance          0.55    0.01 1.00   -1.42   -0.12    0.55    1.23    2.50
## beta_tempo          0.23    0.00 0.04    0.16    0.21    0.23    0.25    0.30
## beta_genre         -0.02    0.00 0.05   -0.12   -0.05   -0.02    0.01    0.07
## lp__             -309.30    0.02 1.88 -313.80 -310.34 -308.97 -307.91 -306.64
##                  n_eff Rhat
## mu_popularity    16326    1
## sigma_popularity 24531    1
## beta_duration    35134    1
## beta_explicito   23850    1
## beta_dance       13663    1
## beta_tempo       21824    1
## beta_genre       22077    1
## lp__             12534    1
## 
## Samples were drawn using NUTS(diag_e) at Thu May 16 00:11:55 2024.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Con base en lo anterior:

La duración no parece tener una influencia significativa en la popularidad, ya que su coeficiente es cercano a cero.
La explicitividad tiene una influencia positiva significativa en la popularidad, con un coeficiente de 1.01, lo que sugiere que las canciones más explícitas tienden a ser más populares.
La danzabilidad también tiene una influencia positiva significativa en la popularidad, con un coeficiente de 0.55, lo que indica que las canciones más adecuadas para bailar tienden a ser más populares.
El tempo también tiene una influencia positiva en la popularidad, aunque más moderada que la explicitividad y la danzabilidad, con un coeficiente de 0.23.
El género parece tener una influencia ligeramente negativa en la popularidad, aunque muy pequeña, con un coeficiente de -0.02, lo que sugiere que ciertos géneros podrían ser ligeramente menos populares en promedio en comparación con otros.

En resumen, la explicitividad, la danzabilidad y el tempo parecen ser los principales impulsores de la popularidad de una canción, mientras que la duración y el género tienen influencias menos significativas.

Observamos los intervalos de credibilidad

extract_lineal <- rstan::extract(fit_lineal)

stan_plot(
  fit_lineal, 
  pars = c("mu_popularity", "beta_duration", "beta_explicito", "beta_dance", "beta_tempo", "beta_genre"), 
  prob = 0.95
)

## ci_level: 0.8 (80% intervals)

## outer_level: 0.95 (95% intervals)

2) Modelo lineal con priors informativas

En nuestro modelo anterior, las prioris eran poco o nulamente informativas por lo que integraremos a dicho modelo prioris informativas que reflejen nuestro conocimiento de las variables:

¿Cuánto dura la canción?: duration_ms
¿Es explicita?: explicit
¿Qué tan bailable es?: danceability
Tempo: tempo
Género: track_genre

Proponemos como distribuciones iniciales las siguientes:

Popularidad:

Es nuestra variable objetivo y debido a que la popularidad es un valor discreto entre 0 y 100, elegimos una distribución Binomial Negativa la cual permitirá modelar la sobredispersión, pues esperamos que la popularidad tenga una amplia variabilidad.

En este caso, consideramos la popularidad como el resultado de un proceso de “éxito”, es decir, la canción es escuchada o marcada como favorita por un usuario hasta que se alcanza un número fijo de “éxitos” (cierto nivel de popularidad).

\[popularidad \sim NegBIn(\mu,\phi)\]
Duración:

Sabemos que la duración de las canciones debe ser mayor a cero y tiene una media aproximada de 3 minutos, por lo que usaremos una distribución exponencial.

\[duración \sim Exp(\lambda) ; \lambda = 1/180000\] (se encuentra en milisegundos)
Explícita:

La variable “explícita” es una variable que tomará los valores 0 y 1. Usaremos entonces una distribución Beta para modelarla.

\[explicita \sim Beta(\rho)\]
Danzabilidad:

Ya que los valores que toma esta variable se encuentran entre 0 y 1, utilizaremos una distribución Beta.

\[danzabilidad \sim Beta(\alpha,\beta)\]
Tempo:

Variable continua entre 0 y 244. Utilizaremos una distribución normal truncada.

\[tempo \sim N+(\mu,\sigma)\]
Género:

El género es una variable entera categórica, por lo que la modelaremos usando una distribución Multinomial.

\[género \sim Multinomial(n,p)\]

Utilizaremos un modelo lineal generalizado (GLM), donde la popularidad sea una función lineal de las demás variables, teniendo en cuenta la naturaleza discreta de la popularidad (valores enteros entre 0 y 100).

  stan_code <- "
  data {
    int<lower=0> N;               
    int<lower=0, upper=100> popularity[N];  
    real<lower=0> duration_ms[N];      
    int<lower=0, upper=1> explicito[N];
    real<lower=0, upper=1> danceability[N];  
    real<lower=0, upper=244> tempo[N];     
    real<lower=1, upper=114> genre[N];        
  }

  parameters {
    real phi_popularity;
    real<lower=0> p_duration;        
    real<lower=0, upper=1> p_explicit;  
    real<lower=0, upper=1> p_dance; 
    real<lower=0> mu_tempo;             
    real<lower=0> sigma_tempo;  
    real<lower=0> mu_genre;             
    real<lower=0> sigma_genre;        
    real<lower=0> p_tempo;                
    real<lower=0> p_genre;              

  }

  model {

    // Priors
    phi_popularity ~ normal(0, 1);
    p_duration ~ exponential(180000);  
    p_explicit ~ beta(2,2);                
    p_dance ~ beta(1, 1);           
    mu_tempo ~ normal(120, 30);               
    sigma_tempo ~ cauchy(0, 5);      
    mu_genre ~ normal(120, 30);               
    sigma_genre ~ cauchy(0, 5);  
    p_tempo ~ normal(mu_tempo, sigma_tempo);
    p_genre ~ normal(mu_genre, sigma_genre);

    // Likelihood
    for (i in 1:N) {
      real mu_popularity;
      
      mu_popularity = p_duration * duration_ms[i] + 
                      p_explicit * explicito[i] +
                      p_dance * danceability[i] +
                      p_tempo * tempo[i] +
                      p_genre * genre[i];
      
      popularity[i] ~ neg_binomial_2(mu_popularity, phi_popularity); // Distribución binomial negativa para popularidad
    }
  }
  "

# Compilar el modelo
stan_model <- stan_model(model_code = stan_code)

stan_data <- list(
  N = nrow(selected_datos_small),
  popularity = selected_datos_small$popularity,
  duration_ms = selected_datos_small$duration_ms,
  explicito = selected_datos_small$explicit,
  danceability = selected_datos_small$danceability,
  tempo = selected_datos_small$tempo,
  genre = selected_datos_small$genre_numeric
)

# Ajustar el modelo a los datos
fit <- sampling(stan_model, data = stan_data, chains = 4, iter = 50000, warmup = 5000, seed=123)

## Warning: There were 31 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.

## Warning: Examine the pairs() plot to diagnose sampling problems

Observamos los resultados

print(fit)

## Inference for Stan model: anon_model.
## 4 chains, each with iter=50000; warmup=5000; thin=1; 
## post-warmup draws per chain=45000, total post-warmup draws=180000.
## 
##                   mean se_mean     sd    2.5%     25%     50%     75%   97.5%
## phi_popularity    0.92    0.00   0.19    0.60    0.79    0.91    1.04    1.33
## p_duration        0.00    0.00   0.00    0.00    0.00    0.00    0.00    0.00
## p_explicit        0.51    0.00   0.22    0.10    0.33    0.51    0.68    0.91
## p_dance           0.50    0.00   0.29    0.03    0.26    0.51    0.75    0.98
## mu_tempo         97.52    0.17  36.55   11.76   75.28   99.64  122.43  163.95
## sigma_tempo     123.27    0.88 260.39   10.98   52.29   82.07  133.61  464.20
## mu_genre         97.19    0.16  36.80   10.31   74.97   99.57  122.25  163.93
## sigma_genre     122.64    0.87 276.36    9.99   51.92   81.45  132.83  460.53
## p_tempo           0.26    0.00   0.06    0.16    0.22    0.26    0.30    0.39
## p_genre           0.12    0.00   0.10    0.00    0.04    0.09    0.16    0.36
## lp__           -264.85    0.02   2.77 -271.41 -266.42 -264.41 -262.82 -260.74
##                 n_eff Rhat
## phi_popularity 149577    1
## p_duration     205032    1
## p_explicit     187765    1
## p_dance        194623    1
## mu_tempo        46470    1
## sigma_tempo     87549    1
## mu_genre        50216    1
## sigma_genre    100424    1
## p_tempo        126896    1
## p_genre        130962    1
## lp__            33605    1
## 
## Samples were drawn using NUTS(diag_e) at Thu May 16 00:12:34 2024.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Podemos deducir que la explicitividad, la danzabilidad, el tempo y el género parecen tener influencias significativas en la popularidad de una canción, mientras que la duración parece tener una influencia insignificante. Específicamente, las canciones más explícitas, adecuadas para bailar y con ciertos tempos y géneros tienden a ser más populares en promedio.

Observamos intervalos de credibilidad

stan_plot(
  fit, 
  pars = c("phi_popularity", "p_duration", "p_explicit", "p_dance", "mu_tempo", "mu_genre", "p_tempo", "p_genre"), 
  prob = 0.95
)

## ci_level: 0.8 (80% intervals)

## outer_level: 0.95 (95% intervals)

3) Modelo de popularidad individual

El siguente modelo considerará solo la popularidad del artista como predictor de la popularidad de la canción.

individual_code <-"
data {
  int<lower=0> N;                
  vector[N] popularity_artist;   
  vector[N] popularity;     
}

parameters {
  real intercept;                 
  real beta_artist;               
  real<lower=0> sigma;            
}

model {
  // Priors
  intercept ~ normal(0, 1);      
  beta_artist ~ normal(0, 1);     
  sigma ~ normal(0, 1);           

  // Likelihood
  for (i in 1:N) {
    // Modelo lineal
    popularity[i] ~ normal(intercept + beta_artist * popularity_artist[i], sigma);
  }
}
"


# Compilar el modelo
individual_model <- stan_model(model_code = individual_code)

stan_indiv_data <- list(
  N = nrow(selected_datos),
  popularity = selected_datos$popularity,
  popularity_artist = selected_datos$popularity_artist
)

# Ajustar el modelo a los datos
individual_fit <- sampling(individual_model, data = stan_indiv_data, chains = 4, iter = 10000, warmup = 500, seed=123)

Observamos los resultados

print(individual_fit)

## Inference for Stan model: anon_model.
## 4 chains, each with iter=10000; warmup=500; thin=1; 
## post-warmup draws per chain=9500, total post-warmup draws=38000.
## 
##                   mean se_mean   sd       2.5%        25%        50%        75%
## intercept         0.00    0.00 0.08      -0.16      -0.06       0.00       0.06
## beta_artist       1.00    0.00 0.00       1.00       1.00       1.00       1.00
## sigma            11.49    0.00 0.03      11.43      11.47      11.49      11.51
## lp__        -220714.03    0.01 1.23 -220717.25 -220714.57 -220713.69 -220713.13
##                  97.5% n_eff Rhat
## intercept         0.17  9286    1
## beta_artist       1.00 11097    1
## sigma            11.54 25286    1
## lp__        -220712.63 12785    1
## 
## Samples were drawn using NUTS(diag_e) at Thu May 16 02:03:28 2024.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

En resumen, el artista sí logra influir en la popularidad de una canción.

Observamos los intervalos de credibilidad

stan_plot(
  individual_fit, 
  pars = c("intercept", "beta_artist"), 
  prob = 0.95
)

## ci_level: 0.8 (80% intervals)

## outer_level: 0.95 (95% intervals)

Para observar el desempeño obtendremos las estimaciones de popularidad para 20 canciones y compararemos estas con la popularidad real.

En el siguiente gráfico puede observarse un intervalo del 95% de credibilidad para las 10 canciones, así como la media de las simulaciones y el valor real de la popularidad para la canción.

avg_popularity_artists <- mean(selected_datos$popularity_artist)

n_songs = 20
# Calculamos las predicciones con los coeficientes obtenidos
comp_simul <- list()

for (i in 1:n_songs) {
  mu_pred <- mean(rstan::extract(individual_fit)$intercept) + 
             mean(rstan::extract(individual_fit)$beta_artist) * avg_popularity_artists
             
  pred_simuladas <- rnorm(10000, mean = mu_pred, sd = mean(rstan::extract(individual_fit)$sigma))

  comp_simul[[i]] <- list(real_popularity = datos$popularity[i], simulated_popularity = pred_simuladas)
}

# Establecer márgenes más pequeños
par(mar = c(2, 2, 2, 2))

par(mfrow=c(5, 4))  

for (i in 1:length(comp_simul)) {

  real_popularity <- comp_simul[[i]]$real_popularity
  simulated_popularity <- comp_simul[[i]]$simulated_popularity
  
  # intervalo de confianza del 95%
  ci <- quantile(simulated_popularity, c(0.025, 0.975))
  

  plot(1, xlim=c(0, 2), ylim=c(0, max(c(real_popularity, ci))), type="n", xlab="", ylab="Popularidad")
  lines(rep(1, 2), ci, col="darkcyan", lwd=2)  
  points(1, abs(mean(simulated_popularity)), col="brown", pch=19)  
  points(1, real_popularity, col="orange", pch=19)  
  text(1.05, real_popularity, paste("Real:", real_popularity), pos=4,offset=1)  # Etiqueta valor real
  text(1.05, mean(simulated_popularity), paste("Media:", round(abs(mean(simulated_popularity)), 2)), pos=4,offset=-6)  # Etiqueta  media
  title(paste("Canción", i))  # Título
}

par(mfrow=c(1, 1)) 

par(mar = c(5, 4, 4, 2) + 0.1)

Observamos que en la mayoría de las canciones que estamos ejemplificando la estimación no está cerca del valor real ni dentro del intervalo de credibilidad.

4) Modelo popularidad entre artistas

general_stan_code <-" 
data {
  int<lower=0> N;                // Número de observaciones
  real popularity_avg_artist;    // Popularidad promedio entre artistas
  vector[N] popularity;          // Popularidad de la canción
}

parameters {
  real intercept;                // Intercepto
  real beta_artist;              // Coeficiente de la popularidad entre artistas
  real<lower=0> sigma;           // Desviación estándar de los errores
}

model {
  // Priors
  intercept ~ normal(0, 1);      // Prior para el intercepto
  beta_artist ~ normal(0, 1);     // Prior para el coeficiente de la popularidad entre artistas
  sigma ~ normal(0, 1);           // Prior para la desviación estándar

  // Likelihood
  for (i in 1:N) {
    // Modelo lineal
    popularity[i] ~ normal(intercept + beta_artist * popularity_avg_artist, sigma);
  }
}

//generated quantities {
//  vector[N] simulated_popularity;  // Simulaciones de la popularidad de las canciones

  // Generar simulaciones
//  for (i in 1:N) {
//    simulated_popularity[i] = normal_rng(intercept + beta_artist * popularity_avg_artist, sigma);
//  }
//}
"


# Compilar el modelo
general_model <- stan_model(model_code = general_stan_code)

avg_popularity_artists <- mean(selected_datos$popularity_artist)

stan_general_data <- list(
  N = nrow(selected_datos),
  popularity = selected_datos$popularity,
  popularity_avg_artist = avg_popularity_artists
)

# Ajustar el modelo a los datos
general_fit <- sampling(general_model, data = stan_general_data, chains = 4, iter = 10000, warmup = 500, seed=123)

Observamos los resultados

print(general_fit)

## Inference for Stan model: anon_model.
## 4 chains, each with iter=10000; warmup=500; thin=1; 
## post-warmup draws per chain=9500, total post-warmup draws=38000.
## 
##                   mean se_mean   sd       2.5%        25%        50%        75%
## intercept         0.02    0.01 1.01      -1.94      -0.65       0.02       0.69
## beta_artist       1.00    0.00 0.03       0.94       0.98       1.00       1.02
## sigma            22.25    0.00 0.06      22.14      22.21      22.25      22.29
## lp__        -270666.10    0.01 1.22 -270669.22 -270666.64 -270665.78 -270665.21
##                  97.5% n_eff Rhat
## intercept         1.99  7023    1
## beta_artist       1.06  7004    1
## sigma            22.36 20275    1
## lp__        -270664.71  9301    1
## 
## Samples were drawn using NUTS(diag_e) at Thu May 16 06:02:46 2024.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Observamos los intervalos de credibilidad

stan_plot(
  general_fit, 
  pars = c("intercept", "beta_artist"), 
  prob = 0.95
)

## ci_level: 0.8 (80% intervals)

## outer_level: 0.95 (95% intervals)

Para observar el desempeño obtendremos las estimaciones de popularidad para 20 canciones y compararemos estas con la popularidad real.

En el siguiente gráfico puede observarse un intervalo del 95% de credibilidad para las 10 canciones así como la media de las simulaciones y el valor real de la popularidad para la canción.

n_songs = 20
comp_simul <- list()

for (i in 1:n_songs) {
  mu_pred <- mean(rstan::extract(general_fit)$intercept) + 
             mean(rstan::extract(general_fit)$beta_artist) * avg_popularity_artists
             
  pred_simuladas <- rnorm(10000, mean = mu_pred, sd = mean(rstan::extract(general_fit)$sigma))

  comp_simul[[i]] <- list(real_popularity = datos$popularity[i], simulated_popularity = pred_simuladas)
}

par(mar = c(2, 2, 2, 2))


par(mfrow=c(5, 4))  

for (i in 1:length(comp_simul)) {

  real_popularity <- comp_simul[[i]]$real_popularity
  simulated_popularity <- comp_simul[[i]]$simulated_popularity
  
  # intervalo de confianza del 95%
  ci <- quantile(simulated_popularity, c(0.025, 0.975))
  

  plot(1, xlim=c(0, 2), ylim=c(0, max(c(real_popularity, ci))), type="n", xlab="", ylab="Popularidad")
  lines(rep(1, 2), ci, col="darkcyan", lwd=2)  
  points(1, abs(mean(simulated_popularity)), col="brown", pch=19)  
  points(1, real_popularity, col="orange", pch=19)  
  text(1.05, real_popularity, paste("Real:", real_popularity), pos=4,offset=1)  # Etiqueta valor real
  text(1.05, mean(simulated_popularity), paste("Media:", round(abs(mean(simulated_popularity)), 2)), pos=4,offset=-6)  # Etiqueta  media
  title(paste("Canción", i))  # Título
}

par(mfrow=c(1, 1)) 

par(mar = c(5, 4, 4, 2) + 0.1)

En este caso, podemos observar que la mayoría de las estimaciones se acercan al valor real o, por lo menos, están dentro del intervalo de credibilidad. Por lo tanto, modelando de esta manera se obtuvo una mejora relativa al modelo anterior.

5) Modelo jerárquico

jerarquico_code <- "
data {
  int<lower=0> N;                     // Número de observaciones
  vector[N] popularity_artist;        // Popularidad individual del artista
  real popularity_avg_artist;         // Popularidad promedio entre artistas
  vector[N] popularity;               // Popularidad de la canción
}

parameters {
  real intercept;                     // Intercepto
  real beta_individual_artist;        // Coeficiente de la popularidad individual del artista
  real beta_avg_artist;               // Coeficiente de la popularidad promedio entre artistas
  real<lower=0> sigma;                // Desviación estándar de los errores
}

model {
  // Priors
  intercept ~ normal(0, 1);          // Prior para el intercepto
  beta_individual_artist ~ normal(0, 1);  // Prior para el coeficiente de la popularidad individual del artista
  beta_avg_artist ~ normal(0, 1);     // Prior para el coeficiente de la popularidad promedio entre artistas
  sigma ~ normal(0, 1);               // Prior para la desviación estándar

  // Likelihood
  for (i in 1:N) {
    // Modelo lineal jerárquico
    popularity[i] ~ normal(intercept + beta_individual_artist * popularity_artist[i] + beta_avg_artist * popularity_avg_artist, sigma);
  }
}
"

# Compilar el modelo
jerarquico_model <- stan_model(model_code = jerarquico_code)

stan_jerarquico_data <- list(
  N = nrow(selected_datos),
  popularity = selected_datos$popularity,
  popularity_avg_artist = avg_popularity_artists,
  popularity_artist = selected_datos$popularity_artist
)

# Ajustar el modelo a los datos
jerarquico_fit <- sampling(jerarquico_model, data = stan_jerarquico_data, chains = 4, iter = 10000, warmup = 500, seed=123)

Observamos los resultados

print(jerarquico_fit)

## Inference for Stan model: anon_model.
## 4 chains, each with iter=10000; warmup=500; thin=1; 
## post-warmup draws per chain=9500, total post-warmup draws=38000.
## 
##                              mean se_mean   sd       2.5%        25%        50%
## intercept                   -0.01    0.01 0.98      -1.92      -0.67      -0.01
## beta_individual_artist       1.00    0.00 0.00       1.00       1.00       1.00
## beta_avg_artist              0.00    0.00 0.03      -0.06      -0.02       0.00
## sigma                       11.49    0.00 0.03      11.43      11.47      11.49
## lp__                   -220714.51    0.01 1.41 -220718.01 -220715.21 -220714.19
##                               75%      97.5% n_eff Rhat
## intercept                    0.65       1.94  6491    1
## beta_individual_artist       1.00       1.00 29721    1
## beta_avg_artist              0.02       0.06  6485    1
## sigma                       11.51      11.55 26403    1
## lp__                   -220713.47 -220712.76 11281    1
## 
## Samples were drawn using NUTS(diag_e) at Thu May 16 09:46:48 2024.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Por lo tanto, la influencia de la popularidad del artista individual parece ser el factor más relevante en la popularidad de una canción, mientras que la influencia del promedio de popularidad de los artistas no parece ser significativa.

Observamos los intervalos de credibilidad

stan_plot(
  jerarquico_fit, 
  pars = c("intercept", "beta_individual_artist", "beta_avg_artist"), 
  prob = 0.95
)

## ci_level: 0.8 (80% intervals)

## outer_level: 0.95 (95% intervals)

Para observar el desempeño obtendremos las estimaciones de popularidad para 20 canciones y compararemos estas con la popularidad real.

En el siguiente gráfico puede observarse un intervalo del 95% de credibilidad para las 10 canciones así como la media de las simulaciones y el valor real de la popularidad para la canción.

n_songs = 20
# Calculamos las predicciones con los coeficientes obtenidos
comp_simul <- list()

for (i in 1:n_songs) {
  mu_pred <- mean(rstan::extract(jerarquico_fit)$intercept) + 
             mean(rstan::extract(jerarquico_fit)$beta_individual_artist) * datos$popularity_artist[i] +
             mean(rstan::extract(jerarquico_fit)$beta_avg_artist) * avg_popularity_artists
             
  pred_simuladas <- rnorm(10000, mean = mu_pred, sd = mean(rstan::extract(jerarquico_fit)$sigma))

  comp_simul[[i]] <- list(real_popularity = datos$popularity[i], simulated_popularity = pred_simuladas)
}

par(mar = c(2, 2, 2, 2))


par(mfrow=c(5, 4))  

for (i in 1:length(comp_simul)) {

  real_popularity <- comp_simul[[i]]$real_popularity
  simulated_popularity <- comp_simul[[i]]$simulated_popularity
  
  # intervalo de confianza del 95%
  ci <- quantile(simulated_popularity, c(0.025, 0.975))
  

  plot(1, xlim=c(0, 2), ylim=c(0, max(c(real_popularity, ci))), type="n", xlab="", ylab="Popularidad")
  lines(rep(1, 2), ci, col="darkcyan", lwd=2)  
  points(1, abs(mean(simulated_popularity)), col="brown", pch=19)  
  points(1, real_popularity, col="orange", pch=19)  
  text(1.05, real_popularity, paste("Real:", real_popularity), pos=4,offset=1)  # Etiqueta valor real
  text(1.05, mean(simulated_popularity), paste("Media:", round(abs(mean(simulated_popularity)), 2)), pos=4,offset=-6)  # Etiqueta  media
  title(paste("Canción", i))  # Título
}

par(mfrow=c(1, 1)) 

par(mar = c(5, 4, 4, 2) + 0.1)

6) Modelo jerárquico incluyendo la valencia con DAG

Incluimos la valencia en los datos

#Variables de interés
full_data <- datos
selected_datos <- datos[, c("popularity", "valence", "duration_ms","explicit","danceability","tempo","genre_numeric","popularity_artist")]
selected_datos$explicit <- as.integer(selected_datos$explicit)
selected_datos <- head(selected_datos,5000)
#selected_datos_small <- head(selected_datos,50)

avg_popularity_artists <- mean(selected_datos$popularity_artist)

Modelo incluyendo la valencia

jer_val_code <- "
data {
  int<lower=0> N;                     // Número de observaciones
  vector[N] popularity_artist;        // Popularidad individual del artista
  real popularity_avg_artist;         // Popularidad promedio entre artistas
  vector[N] popularity;               // Popularidad de la canción
  vector[N] valence;                  // Valencia de la canción
}

parameters {
  real intercept;                     // Intercepto
  real beta_individual_artist;        // Coeficiente de la popularidad individual del artista
  real beta_avg_artist;               // Coeficiente de la popularidad promedio entre artistas
  real beta_valence;                  // Coeficiente de la valencia
  real<lower=0> sigma;                // Desviación estándar de los errores
}

model {
  // Priors
  intercept ~ normal(0, 1);          // Prior para el intercepto
  beta_individual_artist ~ normal(0, 1);  // Prior para el coeficiente de la popularidad individual del artista
  beta_avg_artist ~ normal(0, 1);     // Prior para el coeficiente de la popularidad promedio entre artistas
  beta_valence ~ normal(0, 1);        // Prior para el coeficiente de la valencia
  sigma ~ normal(0, 1);               // Prior para la desviación estándar

  // Likelihood
  for (i in 1:N) {
    // Modelo lineal jerárquico
    popularity[i] ~ normal(intercept + beta_individual_artist * popularity_artist[i] + beta_avg_artist * popularity_avg_artist + beta_valence * valence[i], sigma);
  }
}
"

jer_val_model <- stan_model(model_code = jer_val_code)

stan_jer_val_data <- list(
  N = nrow(selected_datos),
  popularity = selected_datos$popularity,
  valence = selected_datos$valence,
  popularity_avg_artist = avg_popularity_artists,
  popularity_artist = selected_datos$popularity_artist
)

# Ajustar el modelo a los datos
jer_val_fit <- sampling(jer_val_model, data = stan_jer_val_data, chains = 4, iter = 10000, warmup = 500)

Observamos los resultados

print(jer_val_fit)

## Inference for Stan model: anon_model.
## 4 chains, each with iter=10000; warmup=500; thin=1; 
## post-warmup draws per chain=9500, total post-warmup draws=38000.
## 
##                             mean se_mean   sd      2.5%       25%       50%
## intercept                   0.00    0.01 0.99     -1.95     -0.66      0.00
## beta_individual_artist      1.00    0.00 0.01      0.99      1.00      1.00
## beta_avg_artist             0.00    0.00 0.03     -0.06     -0.02      0.00
## beta_valence                0.07    0.00 0.40     -0.72     -0.20      0.06
## sigma                       8.00    0.00 0.08      7.85      7.95      8.00
## lp__                   -12960.92    0.01 1.57 -12964.81 -12961.73 -12960.60
##                              75%     97.5% n_eff Rhat
## intercept                   0.66      1.94 16785    1
## beta_individual_artist      1.00      1.01 36961    1
## beta_avg_artist             0.02      0.06 16415    1
## beta_valence                0.34      0.86 22370    1
## sigma                       8.06      8.16 27633    1
## lp__                   -12959.77 -12958.85 15587    1
## 
## Samples were drawn using NUTS(diag_e) at Thu May 16 11:01:22 2024.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

La influencia de la popularidad del artista individual parece ser el factor más relevante en la popularidad de una canción, mientras que la influencia del promedio de popularidad de los artistas y la valencia también son significativas, pero en menor medida.

Observamos los intervalos de credibilidad

stan_plot(
  jer_val_fit, 
  pars = c("intercept", "beta_individual_artist", "beta_avg_artist", "beta_valence"), 
  prob = 0.95
)

## ci_level: 0.8 (80% intervals)

## outer_level: 0.95 (95% intervals)

DAG (Grafo Acíclico Dirigido, por sus siglas en inglés)

*Popularidad de la canción (popularity): La popularidad de una canción puede ser influenciada por diversos factores, entre ellos, la calidad de la música en sí misma y cómo resuena con el público objetivo.

*Popularidad individual del artista (popularity_artist): La popularidad de un artista puede influir en la percepción y el interés del público hacia sus nuevas canciones. Los seguidores de un artista pueden estar más dispuestos a escuchar y apoyar nuevas producciones.

*Popularidad promedio entre artistas (popularity_avg_artist): La popularidad promedio entre artistas podría representar el contexto musical general en el que se encuentra la canción. Por ejemplo, si la mayoría de las canciones lanzadas recientemente por otros artistas tienen una alta popularidad, esto podría aumentar las expectativas del público y afectar la percepción de la nueva canción.

*Valencia (valence): La valencia de la música, es decir, su capacidad para transmitir emociones y sensaciones auténticas, puede tener un efecto directo en cómo el público la percibe y la valora. Una música con alta valencia podría resonar más con los oyentes y aumentar su disposición a compartir y apoyar la canción.

Proponemos un DAG donde las flechas indican la dirección de la causalidad, sugiriendo que la popularidad del artista y la popularidad promedio entre artistas podrían influir en la valencia de la música, que a su vez afecta la popularidad de la canción.

grViz("
digraph {
  graph [ranksep = 0.2, rankdir = LR]
  node [shape=plaintext]
  C
  A
  P
  V
  edge [minlen = 3]
  A -> P -> C
  V -> C
}
", width = 400)

donde C es la popularidad de la canción, A es la popularidad individual del artista, P es la popularidad promedio entre artistas y V es la valencia de la canción.

Resultados

Análisis general de popularidad mediante Modelos Lineales Generalizados

Los tres modelos iniciales creados para analizar la popularidad de canciones en Spotify ofrecen un panorama progresivo de cómo diversos factores influyen en la percepción y aceptación de una canción por parte de los usuarios de la aplicación. El primer modelo, basado en priors no informativas, establece un punto de partida fundamental al recopilar cinco variables predictoras clave; este enfoque inicial permite comprender la relación entre estas variables y la popularidad de las canciones. En contraste, el segundo modelo introduce priors informativas que aprovechan el conocimiento previo para mejorar la precisión de las predicciones. Esta inclusión de información contextual refleja de mejor forma la realidad y las tendencias de popularidad de las canciones. Finalmente, el tercer modelo expande aún más la comprensión al agregar una variable adicional, lo que permite evaluar cómo este nuevo factor afecta la popularidad de las canciones. Este enfoque más completo revela cómo la inclusión de diferentes aspectos puede modular la percepción de una canción por parte del público.

En resumen, estos modelos proporcionan una perspectiva incremental y enriquecedora para entender la dinámica detrás de la popularidad de las canciones en la plataforma de Spotify, destacando la importancia de considerar tanto los datos como el contexto para obtener resultados precisos.

Modelos heterogéneo, homogéneo y jerárquico

Los modelos heterogéneo y homogéneo planteados en este trabajo, proporcionan perspectivas simples pero esclarecedoras. En el primer modelo calculamos la media de la popularidad para cada artista de manera individual, lo que ofreció una visión detallada de cómo cada artista impacta la popularidad de sus canciones sin considerar la interacción con otros artistas. El segundo modelo, al calcular la media de todos los artistas de forma conjunta, ofrece una visión generalizada al considerar la interacción de la popularidad entre todos los artistas, sin embargo, perdimos la distinción de la influencia de cada artista en particular.

En contraste, el modelo jerárquico adopta una perspectiva más compleja, pues considera tanto la variabilidad entre artistas como la media general de popularidad. Este enfoque nos permitió capturar tanto la influencia individual de cada artista como la variabilidad global entre ellos, ofreciéndonos una comprensión más completa del fenómeno.

Finalmente, el modelo jerárquico incluyendo la valencia nos permitió considerar esta variable cómo explicativa de la popularidad de una canción. Esto es relevante ya que consideramos, como se puede ver en el DAG, que la valencia influye en la popularidad de una canción. Este efecto se suma al de la popularidad individual del artista que a su vez impacta la popularidad promedio y determinan la popularidad de una canción.

Así, aunque los modelos planos simplifican el análisis al centrarse exclusivamente en el efecto del artista, el enfoque jerárquico proporcionó una visión más completa al permitir la evaluación de múltiples niveles de influencia. Sin embargo, únicamente incluimos la popularidad dada, las popularidades generadas individuales y colectivas de los artistas y la valencia, por lo que consideramos que la exclusión de otras variables predictoras pudo haber limitado la comprensión total de la popularidad de las canciones en Spotify. Por lo tanto, en futuros análisis consideraremos la inclusión de más variables predictoras para obtener una imagen más precisa y completa del problema.

Conclusiones

En conclusión, aunque los modelos planos simplifican el análisis al centrarse exclusivamente en el efecto del artista, el enfoque jerárquico proporcionó una visión más completa al permitir la evaluación de múltiples niveles de influencia. Sin embargo, únicamente incluimos la popularidad dada, las popularidades generadas individuales y colectivas de los artistas y la valencia, por lo que consideramos que la exclusión de otras variables predictoras pudo haber limitado la comprensión total de la popularidad de las canciones en Spotify. Por lo tanto, en futuros análisis consideraremos la inclusión de más variables predictoras para obtener una imagen más precisa y completa del problema.

Fuentes

Johnson, A. A., Ott, M. Q., & Dogucu, M. (2022). Bayes rules! an introduction to applied Bayesian modeling. CRC Press.
Kaggle. (2021). Spotify Tracks Dataset [Archivo de datos CSV]. Recuperado de https://www.kaggle.com/datasets/maharshipandya/-spotify-tracks-dataset?select=dataset.csv

Anexos

Anexo 1: Limpieza, Adecuación y Exploración de la Base de Datos de Spotify

El Anexo 1 incluye información detallada sobre la exploración inicial de la base de datos de canciones de Spotify, así como los procedimientos de limpieza y adecuación de los datos. Se presenta un análisis exhaustivo de la estructura y calidad de los datos, identificando y abordando posibles errores, valores atípicos y datos faltantes. Además, se describen las transformaciones realizadas para preparar los datos para su posterior análisis.

Anexo 2: Modelos en cmdstanr

El Anexo 2 incluye otros modelos donde se compara la popularidad promedio de canciones para un mismo artista.